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BACKGROUND OF THE INVENTION 

The present invention is related to the field of molecular modeling and, 
more particularly, to computer-implemented methods for the dynamic modeling and 
static analysis of large molecules. 

The field of molecular modeling has successfully simulated the motion 
(molecular dynamics or (MD)) and determined energy minima or rest states (static 
analysis) of many complex molecular systems by computers. Typical molecular 
modeling applications have included enzyme-ligand docking, molecular diffusion, 
reaction pathways, phase transitions, and protein folding studies. Researchers in the 
biological sciences and the pharmaceutical, polymer, and chemical industries are 
beginning to use these techniques to understand the nature of chemical processes in 
complex molecules and to design new drugs and materials accordingly. Naturally, the 
acceptance of these tools is based on several factors, including the accuracy of the 
results in representing reality, the size and complexity of the molecular systems that 
can be modeled, and the speed by which the solutions are obtained. Accuracy of 
many computations has been compared to experiment and generally found to be 
adequate within specified bounds. However, the use of these tools in the prior art has 
required enormous computing power to model molecules or molecular systems of 
even modest size to obtain molecular time histories of sufficient length to be useful. 

There are two sources of computational complexity for molecular 
modeling tasks involving time integration: 

1 . The particular molecular model which is used to describe the locations, 

velocities and mass properties of the constituent atoms, the inter-atomic forces 



between them, and the interactions between the atoms and their surrounding 
environment; and 

2. The particular numerical method used to advance the model through time. 
Time is advanced repeatedly by very short intervals, called timesteps, until a 
5 final time has been reached. 

Substantial work has been completed in reducing the computational 
load for molecular models, such as the reduction of model complexity by constraining 
higher order modes with rigid body assumptions, simplifying the model with rigid or 
flexible substructuring, Order(7V) dynamics, efficient implicit solvent models, and 
10 multipole methods for the force field models (see, for example, U.S. Patent No. 
5,424,963 on the commercial MBO(N)D software package). Co-pending 

applications, U.S. Appln. No. , entitled "METHOD FOR LARGE 

C TIMESTEPS IN MOLECULAR MODELING," and U.S. Appln. No. , 

~ Z entitled "METHOD FOR RESIDUAL FORM IN MOLECULAR MODELING," both 

tf'l 

yj 15 of which are filed of even date, claim priority from the previously cited provisional 
*i patent applications and which are assigned to the present assignee, and which are 

incorporated by reference herein, describe further improvements in molecular models 
and numerical methods. 

The primary reason molecular simulations are so slow is that current 
20 numerical methods require very small timesteps, typically between 1 and 10 

femtoseconds (10" 15 to 10" 14 seconds). Each timestep taken requires the computation 
of a new state (position and motion for each atom) of the particular molecular model, 
and then computation of the new set of forces resulting from the new state. For 
example, molecular dynamics simulations of the complex behavior of large 
25 molecules, such as the folding of a protein, typically need to cover a time span from at 
least a microsecond up to several seconds or even minutes. With techniques currently 
in common use, this results in the requirement to take 10 9 to 10 16 timesteps in the 
computer simulation. The per-step computations of the state, and especially the 
forces, grow very expensive as the problem size increases. Even with the fastest 
30 computers available today, months, years or even centuries of computer time are 
required to solve such problems even for systems of modest size. 

Heretofore, it has been widely believed by molecular dynamicists that 
these small timesteps are an inevitable requirement of the need to maintain accuracy 
in the presence of the very high frequencies to be found in vibrations of molecular 
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bonds. For example, see Leach, Molecular Modelling Principles and Applications, 
1996, p. 328; Berendsen, in Computational Molecular Dynamics: Challenges, 
Methods, Ideas Deuflhard et al. (ed.), Springer, 1999, pg. 18; Rapaport, The Art of 
Molecular Dynamics Simulation, Cambridge, 1995, reprinted with corrections 1998, 
5 p. 57; and U.S. Patent No. 5,424,963. 

This common-sense belief is incorrect, however. The computer 
science sub-discipline of numerical analysis has produced an extensive theory of 
numerical integration for problems in which high frequencies exist but are of little 
interest. These problems are termed "stiff problems (see, for example, Hairer and 

1 0 Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic 
Problems, 2 nd ed., Springer, 1996). In these cases, it is the stability of the integration 
method, not the required solution accuracy, which limits the timestep. Integration 
methods exist which have unconditional stability, meaning that in theory they can 
take arbitrarily large timesteps. Such methods have a mathematical property called 

1 5 "L-stability." Only so-called "implicit" integration methods can be L-stable, but very 
few implicit integration methods actually are L-stable. Previously cited co-pending 

U.S. Patent Appln. No. , entitled "METHOD FOR LARGE TIMESTEPS 

IN MOLECULAR MODELING," covers the use of implicit and in particular L-stable 
integration methods. 

20 The present invention covers another critical aspect of allowing the 

implicit and in particular L-stable integrators to take large timesteps, namely, more 
accurate methods for computing required components of the implicit integration 
methods called "Jacobians". 

But the same problem of high frequency molecular vibration for MD 

25 simulations causes problems for the calculation of Jacobians. For example, the 
repulsive van der Waal's forces that are generated as the electron orbitals of two 
atoms begin to overlap must be accounted for in a molecule. This quantum 
mechanical effect is often approximated by the so-called Lennard- Jones potential 
(Rapaport, op. cit.), which models the repulsive force as being proportional to 1/r 13 , 

30 where r is the distance between the centers of the atoms. This is an extreme stiff 

interatomic force which is characteristic of molecular dynamics (MD) simulations and 
poses particular difficulty for any numerical integration scheme used to advance time 
during the simulation. As a result and as mentioned previously, most prior art MD 
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integration schemes have timesteps limited by the high frequencies generated by these 
extremely stiff repulsive forces. And in particular, the stiffness of the atomic forces 
greatly magnify the difficulty of forming certain required Jacobian matrices. Such 
Jacobian matrices are a necessary ingredient of any stable implicit integration scheme, 
such as described in the immediately cited co-pending application. 

All general-purpose simulation codes provide routines to compute 
Jacobians numerically as follows. For a given equation to integrate y = f(y, i) , the 
desired Jacobian is J = df Idy and is numerically computed: 

Ay 

where 

a/ = /U-/U, 

Note that the perturbation Ay has to be selected to give a good result 
and may be difficult to choose, especially for stiff systems. In contrast, analytic 
Jacobians are computed by solving directly, or in this case algorithmically, for the 
equation of the desired derivatives. 

Linearized models are regularly produced analytically for simple 
systems. Such linearization is usually performed around an equilibrium solution. It is 
common in such packages as ACSL (Advanced Continuous Simulation Language), 
(ACSL Reference Guide, Mitchell Gauthier and Associates, 1996), and SPICE (a 
circuit simulation package), (R. Kielkowski, Inside SPICE, McGraw-Hill, 1998) to 
perform equilibrium analysis followed by linearization. Such linearization is 
performed numerically. 

In contrast, the Jacobian of the present invention represents 
linearization about an instantaneous solution of the differential equations (non- 
equilibrium) and is generated analytically. It should be noted that another prior art 
approach to generating analytic Jacobians is to use automated differentiation tools 
such as ADIFOR (Automatic Differentiation of Fortran) (C. Bischof, et. al., ADIFOR 
2.0 Users' Guide, Argonne National Laboratory, 1998) that can symbolically 
differentiate arbitrary equations. These tools could be used to implement this 
invention in practice. However, the structure of the equations must be exploited 
properly to minimize the computations, to avoid errors due to round off and term 
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cancellations, and to avoid "equation swell" which could limit the size of problems 
solved. 

Rather, the present invention allows for the calculation of analytic 
Jacobians for the effective implicit integration, including L-stable integrators, of the 
5 equations of motion of molecular models. 



O 



SUMMARY OF THE INVENTION 
The present invention provides for a method of modeling the behavior 
of a molecule. The method has the steps of selecting a torsion angle, rigid multibody 
10 model for said molecule, the model having equations of motion; selecting an implicit 
integrator; and generating an analytic Jacobian for the implicit integrator to integrate 
the equations of motion so as to obtain calculations of the behavior of the molecule. 
The analytic Jacobian J is derived from an analytic Jacobian of the Residual Form of 
the equations of motion and is described as: 



dq 


du 


du 


du 


dq 


du 



" 9 dq dq m du du 

where q are the generalized coordinates, u are the generalized speeds, 
W is a joint map matrix and M is the mass matrix and p u is the dynamic residual of the 
equations of motion, and z is -M' 1 p u (q,u,0) . The method can also be used for any 
physical system which can be modeled by a torsion angle, rigid multibody system. 

The present invention also provides for the computer code for 
simulating the behavior of a molecule, or any physical system, which can be modeled 
by a torsion angle, rigid multibody system. A module in the computer code with an 
implicit integrator includes the analytic Jacobian as described above. 
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BRIEF DESCRIPTION OF THE DRAWINGS 
Fig. 1 is a representational block module diagram of the software 
system architecture in accordance with the present invention; 

Fig. 2 illustrates the tree structure of the multibody system of the 
molecular model according to the present invention; 

Fig. 3 illustrates the reference configuration of the Fig. 2 multibody 

system; 

Fig. 4A illustrate a sliding joint between two bodies of the Fig. 2 
multibody system; Fig. 4B illustrate a pin joint between two bodies of the Fig. 2 
multibody system; Fig. 4C illustrate a ball joint between two bodies of the Fig. 2 
multibody system; 

Fig. 5 summarizes general computational steps for the Residual Form 
method and Direct Form methods used for the analytic Jacobian computation; 

Fig. 6 is a chart which summarizes the general computational steps for 
7.1 15 the analytic Jacobian; 

Fig. 7 is a plot of digits of accuracy versus perturbation to show the 
accuracy of analytic Jacobian over the numerical Jacobian. 

DESCRIPTION OF THE SPECIFIC EMBODIMENTS 

£ 

RJ 20 GENERAL DESCRIPTION OF THE PRESENT INVENTION 

o 

fc* The numerical method used to advance time in the simulation of a 

modeled molecular system is called the integrator, or integration method. The 
integration proceeds by discretizing the governing equations of motion of the 
molecular system. In the case of an implicit integrator, this step results in a set of 

25 coupled nonlinear algebraic equations (the "stage" equations) and the solution of 

these equations yields the state vector for the molecular system at the next time step. 

The present invention aids the solution of the stage equations. Because 
the atomic forces vary so rapidly over short distances, it is difficult to propagate 
atomic trajectories accurately. Small errors in the atomic positions lead to gross 

30 errors in the satisfaction of the stage equations. Since the Jacobian is used solve the 
stage equations iteratively, an inaccurate Jacobian leads to trial solution that are far 
from the correct solution. This leads to retraction of trial solutions and hinders the 
simulation. 
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Numerical Jacobians may be correct in only half their significant 
digits. An analytical Jacobian will often be correct in all but the last digit. The 
benefit of this result is that the integrator can take far fewer time steps to simulate the 
specified interval, allowing full exploitation of the theoretical stability of the 
integration method. 

The ease or difficulty in producing the Jacobian depends crucially 
upon the formulation used to produce the governing equations. For instance, global 
Cartesian formulations produce equations with very limited explicit coordinate 
dependence. Producing an analytic Jacobian for such a formulation is well known. 
On the other hand, using internal coordinates (in which each molecular subunit's 
location is expressed relative to an earlier subunit's location) as independent variables 
has great benefits. This method is most valuable for any molecule modeled with any 
choice of internal coordinates, and in particular when used with protein models or 
other polymeric molecules using "torsion angles" between the residues as internal 
coordinates. An algorithm for producing an analytic Jacobian for a system formulated 
in torsion angles is extremely difficult to devise. However, the present invention 
achieves this task. 

The present invention addresses a seemingly intractable problem: 
production of the analytic Jacobian for a formulation using internal coordinates, and 
fij 20 specifically torsion angles, which is generally thought to be impractical. In addition 
to being more accurate than numerical Jacobians, the analytic Jacobians are also 
cheaper (in computing power) to produce. The present invention also recognizes a 
key result that the Jacobian of the state derivatives can be computed by applying a 
matrix inverse to the Jacobian of the computed torque method. This result allows 
25 significant savings in computer time and effort to construct the algorithm. 

Furthermore, this method of producing analytic Jacobians for 
multibody system formulations using torsion angle, internal coordinates has not been 
seen in the general MBS literature. The present invention can be used for any torsion 
angle MBS formulation, which can be applied to many other disciplines besides 
30 molecular simulations, including, but not necessarily limited to, mechanical systems, 
robotic systems, vehicle systems, or any other system which could be described as a 
set of hinge-connected rigid bodies. 

OVERVIEW OF DESCRIPTION 
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The preferred embodiment is divided into several sections. The first 
set of sections describes the MD simulation architecture, the multibody system (MBS) 
definitions, and Residual Form of the MBS equations for the subsequent descriptions. 
The next set of sections discusses the definition of the Jacobian, its role in the implicit 
5 integration method, and the computation of the analytic Jacobian using the Residual 
Form. Also shown is the superior accuracy and performance of the analytic Jacobian 
vs. the numerical Jacobian. Further efficiencies in the computation of the analytic 
Jacobian are discussed, specifically, exploiting the rigid body MBS to "contract" the 
size of the Jacobian matrix, and exploiting the topological structure of the MBS to 
10 eliminate unnecessary computations. 

To solve ordinary differential equations (ODEs), most of the prior art 
have used equations expressed in the Direct Form, i.e., y = f {y,t) (where y means 

~~)- The equations of motion for a biomolecular system can be cast into this form 

(and called the Direct Form). In molecular modeling, all prior art known to the 
1 5 present inventors have used the Direct Form. That is, q-Wu, u = M~ x f , where q 
and u are generalized coordinates and speeds respectively, so that conventional ODE 
solution methods can be applied. However, this requires a matrix inversion of 
M (representing the mass of the system) at a cost of Order( N ) to Order( iV 3 ) floating 
point operations (depending on algorithm used, where N is the number of degrees of 
20 freedom in the system), since the natural form of the equations gives rise to inertial 
coupling between the derivatives of the generalized speeds. That is, the equations are 
most naturally produced in the form q-Wu = 0, Mil -f = 0, where M , the mass 
matrix, depends explicitly upon the generalized coordinates q, i.e., M = M(q) . This 
fact requires forming and effectively factoring the mass matrix each time the state 
25 derivatives are needed by the integrator in integrating the equations of motion over 
time. The generalized joint map matrix W is block diagonal and, although it is also 
dependent on the coordinates W = W (q) , it does not have a significant computational 
cost. 

In accordance with the present invention, a method for the solution of 
30 the equations of a molecular system is expressed in Residual Form to bypass the 

customary step of producing the state derivatives directly. The Residual Form method 
has the following steps: 
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1) Discretization of the solution variables. The specific form of discretization is 
dictated by the particular implicit integration method used to advance the 
molecular model in time. Implicit integration follows from the Residual 
Form. Implicit integration, especially L-stable integrators and other highly 
stable integrators, such as implicit Euler, Radau5, SDIRK3, SDIRK4, other 
implicit Runge-Kutta methods, and DASSL or other implicit multistep 
methods, also provide other advantages for molecular modeling. See, for 

example, the above-cited U.S. Patent Appln. No. , entitled "METHOD 

FOR LARGE TIMESTEPS IN MOLECULAR MODELING," filed of even 
date. As a particularly simple example, when used with implicit Euler 
integration, the discretization is as follows: 

q = (q n -q„-i)!h, u = (u n -u n _ 1 )/h 

where h is the timestep. 

2) Substitution into the residual equations: 



pj [M(q)u-f(t,q,u)) 

ft>°< 



s 3) Solution of the resulting nonlinear algebraic equations q = 0 for q n and 

M= 
fU 

ro 

The kinematic residual p q compares an estimated q generated from 

I* 

the implicit integrator to the derivatives computed by the routines for determining the 
20 joints of the molecular model, which is described in greater detail below. The second 

row of the residual is p u , the dynamic residual, which determines the degree to which 

an estimated u satisfies the equations of motion. 

The system mass matrix M and the so-called 'bias-free hinge torque' 

/ are both state dependent. The bias-free hinge torque is generated by the dynamic 
25 residual routine when the calculated u vector passed to the residual routine is zero. In 

general, the hinge accelerations are a response to applied forces, joint torques, and 

motion-induced effects (such as Coriolis and centrifugal forces.) If the system were at 

rest, and subjected only to joint torques, it would be considered in a bias-free state. 

The real system with its actual inputs can be reduced to a bias-free state by computing 
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mathematically. At the core of the module 54 is a multibody system submodule 66. 
The analysis module 56, which communicates with the physical model module 54 and 
the visualization module 58, provides solutions to the computational models of the 
molecular systems defined by the physical model module 54. The analysis module 56 
5 consists of a set of integrator submodules 68 which integrate the differential equations 
of the physical model module 54. The integrator submodules 68 advance the 
molecular system through time and also provide for static analyses used in 
determining the minimum energy configuration of the molecular system. The 
analysis module 56 and its integrator submodules 68 contain most of the subject 
10 matter of the present invention and are described in detail below. 

The visualization module 58 receives input information from the 
biochem components module 52 and the analysis module 56 to provide the user with a 
*f three-dimensional graphical representation of the molecular system and the solutions 

c 

61 obtained for the molecular system. Many visualization modules are presently 

sjj h 15 available, an example being VMD (A. Dalke, et ah, VMD User's Guide . Version 1.5, 

iU June 2000, Theoretical Biophysics Group, University of Illinois, Urbana, Illinois). 

gj The described software code is run on conventional personal 

f computers, such as PCs with Pentium III or Pentium IV microprocessors 

H> manufactured by Intel Corporation of Santa Clara, California. This contrasts with 

|jj 20 many current efforts in molecular modeling which use supercomputers to perform 

calculations. Of course, further speed improvements can be obtained by running the 

described software on faster computers. 

MOLECULAR MODEL AND MULTIBODY SYSTEM DESCRIPTION 

The integrators described below in the submodule 68 operate upon a 

25 set of equations which describe the motion of the molecular model in terms of a 
multibody system (MBS). To aid the computation of the integration methods 
described in detail below, a torsion angle, rigid body model is used to describe the 
subject molecule system, in accordance with the present invention. Internal 
coordinates (selected generalized coordinates and speeds) are used to describe the 

30 states of the molecule. 

The MBS is an abstraction of the atoms and effectively rigid bonds that 
make up the molecular system being modeled and is selected to simplify the actual 
physical system, the molecule in its environment, without losing the features 
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important to the problem being addressed by the simulation. With respect to the 
general system architecture illustrated in Fig. 1, the MBS does not include the 
electrostatic charge or other energetic interactions between atoms nor the model of the 
solvent in which the molecules are immersed. The force fields are modeled in the 
5 submodule 62 and the solvent in the submodule 64 in the biochem components 
module 52. 

Fig. 2 illustrates the tree structure of the MBS of a subject molecule. 
The basic abstraction of the MBS is that of one or more collections of hinge- 
connected rigid bodies 170. A rigid body is a mathematical abstraction of a physical 
1 0 body in which all the particles making up the body have fixed positions relative to 
each other. No flexing or other relative motion is allowed. A hinge connection is a 
mathematical abstraction that defines the allowable relative motion between two rigid 
~I bodies. Examples of these rigid bodies and hinge connections are described below. 

Q One or more of the bodies, called base bodies 172, have special status 

Jy 15 in that their kinematics are referenced directly to a reference point on ground 1 74. 
'% The system graph is one or more "trees". An important property of a tree is that the 

X path from any body to any other body is unique, i.e., the graph contains no loops. The 

Ik bodies in the tree are n in number (the base has the label 1). The bodies in the tree are 

Jf assigned a regular labeling, which means that the body labels never decrease on any 

III 20 path from the base body to any leaf body 1 76. A leaf body is one that is connected to 
2 only a single other body. A regular labeling can be achieved by assigning the label n 

to one of the leaf bodies 178 (there must be at least one). If this body is removed 
from the graph, the tree now has n-l bodies. The label n - 1 is then assigned to one 
of its leaf bodies 180, and the process is repeated until all the bodies have been 
25 labeled. This is also done for any remaining trees in the system. 

To help maintain the relationship between the bodies, an integer 
function is used to record the inboard body for each body of the system. The inboard 
body for each base is ground and i, the parent or inboard body 182 for body k 184, is 
referred to as i = inb(k) . Additionally, the symbol N refers to the inertial, or ground 
30 frame 174. A superscript O refers to the ground origin (0,0,0). 

The symbol for the vector from one point to another contains the name 
of the two points. Thus, r PQ is the vector from the point P to point Q. A vector 
representing the velocity of a point in a reference frame contains the name of the point 
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and the reference frame: N v p . Certain symbols to be introduced later relate two 
reference frames. In this case, the symbol contains the name of two frames. Thus, 
'C k is the direction cosine matrix for the orientation of frame k in frame i. This 
symbol refers to the direction cosine matrix for a typical body in its parent frame. 
5 Thus, l C k (j) indicates the actual body j in question. The left and right superscripts do 
not change with the body index. This is also true for the other symbols. An asterisk 
indicates the transpose: H*(k) , for example. A tilde over a vector indicates a 3 by 3 
skew-symmetric cross product matrix: vw = v x w . £ is an i by i identity matrix., and 
is a zero vector of length i and £■ is an i by i zero matrix. 

1 0 Rigid Bodies of the Model 

Fig. 3 illustrates the reference configuration 190 of a sample "tree" of 
the MBS. More than one tree is allowed. A point of each body is designated as Q, its 
hinge point. For example point Q k 186 is the hinge point for body & 184. A fixed set 
of coordinate axes is established in the inertial frame 198. An arbitrary configuration 

15 of the MBS is chosen as its reference configuration 190. While in this configuration 
the image of the inertial coordinate axes is used to establish a set of body-fixed axes 
in each body. In the reference configuration each hinge point Q is coincident with P, 
a point of its parent body (or extended body.) For each body, point P is called the 
body's inboard hinge point. So the inboard hinge point P k 188 for body k 184 is a 

20 point fixed in its parent body i 182. The inboard hinge point for each base body is a 
point O 192 fixed in ground. The expanded view that shown in Fig. 2 more clearly 
shows that point Q k 1 86 is fixed in body £ 1 84 and point P k 1 88 is fixed in parent 
body / 182. 

The hinge point locations define d(k) 194, a constant vector for each 
25 body, and can also be written r Q ' p " . The vector for body k is fixed in its parent body i. 
It spans from the hinge point for body i to the inboard hinge point for body k. The 
vector d(l) 196 spans from the inertial origin to the first base body's inboard hinge 
point (also a point fixed in ground), and can be written r OQl . 

For a body, m(k) , p(k) , and }g k (k) define the mass properties of 
30 body k for its hinge point Q k . These are, respectively, the mass, the first mass 

moment, and the inertia matrix of the body for its hinge point in the coordinate frame 
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of the body. For a rigid body made up of a distribution of particles, the mass 
properties are constants that are computed by a preprocessing module. The details of 
these computations can be found in standard references, such as Kane, T.R., 
Dynamics, 3 rd Ed., January 1978, Stanford University, Stanford, CA. 
5 Let M(k) , the spatial inertia of body k for its hinge point Q k , be given 

by the symmetric 6 by 6 matrix 

Each joint in the system is described by geometric data. For instance, 
a pin joint is characterized by an axis fixed in the two bodies connected by the joint. 
1 0 The particular data for a joint depends on its type. The number n, the inb function, the 
system mass properties, the vectors d(k) , and the joint geometric data (including joint 
type) constitute the system parameters. 

Joints and Generalized Coordinates of the Model 

Figs. 4A-4C illustrate the joint definitions of the preferred embodiment 

15 of the MBS: the slider joint 100, the pin joint 102, and the ball joint 104. Each joint 
allows translational or rotational displacement of the hinge point Q k 106 relative to 
the inboard hinge point P k 108. These displacements are parameterized by q(k) 110, 
the generalized coordinates for body k. In passing, it should be noted that generalized 
coordinates are examples of generalized quantities, which refer to quantities that have 

20 both rotational character and translational character. For instance, a generalized force 
acting at a point consists of both a force vector and a torque vector. The generalized 
coordinate q(k) for the slider joint 100 is the sliding displacement x 112. The 
generalized coordinate q(k) for the pin joint 102 is the angular displacement 9 1 14. 
The generalized coordinate q(k) for the ball joint 104 is the Euler parameters 

25 (s l ,s 2 ,s 3 ,s 4 ) 116. 

Each joint may be a pin, slider, or ball joint; or a combination of these 
joints. Many other joint types are possible, including, but not limited to, free joints, 
U-joints, cylindrical joints, and bearing joints. For instance, q(k) = (x,y,z), the 
inertial measure numbers of the vector from the base body inboard hinge point to the 

30 base body hinge point express the base body displacement in ground as three 
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orthogonal slider joints. A free joint consists of three orthogonal slider joints 
combined with a ball joint, and has the full 6 degrees of freedom. 

The collection of generalized coordinates for all the bodies comprises 
the vector q , the generalized coordinates for the system. 

Given the generalized coordinates for a particular joint, two quantities: 
r PkQk (k) , the joint translation vector and l C k (k) , the direction cosine matrix for body 
k in its parent are formed. The translation vector r PtQk (k) expresses the vector from 
the inboard hinge point P of body k to the hinge point Q of body k, in the coordinate 
frame of the parent body. Details of these computations depend on the joint type and 
can be easily derived. For purposes of this description, access to a function that can 
generate r PkQk (k) and 'C k (k) given the system generalized coordinates is assumed. 



However, judicious choice greatly simplifies matters. For instance, for pin joints the 
hinge point should be chosen as a point on the axis of the joint. For this choice points 
P and Q remain coincident for all values of the joint angle, so the joint translation is 
zero. If the point Q is chosen at a distance from the axis, points P and Q move 
relative to each other: 



where kis the joint axis unit vector, 0is the joint angle, and r OQk is the vector from 
any point on the axis to point Q. 

For pin joints and ball joints, a point on the axis is always chosen as 
the hinge point. For these joints the translation vector r PkQk (k) is zero. 

For a slider joint, the translation vector r PkQk (k) is q(k)'k . 

The direction cosine matrix for a pin is 

'C* (k) = E 3 cos0+i sin 9 + %X* (1 - cos 0) . 
The direction cosine matrix for a slider is E_ 3 . 

Generalized Speeds of the Model 

Let l V k {k) , the generalized velocity of the hinge point of body k 
measured in its parent i, be parameterized by u(k) , a set of generalized speeds. Then: 



As introduced, the choice of hinge point for each body is arbitrary. 



r PkQk (k) = kx r OQk sin 0-{\- cos 0) (E 3 - U,* ) r c 



OQ k 
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Here, the matrix H(k) is called the joint map for this joint. It is a n u (k) by 6 matrix, 
where n u (k) is the number of degrees of freedom for the joint (1 for a pin or slider, 3 
for a ball, 6 for a free joint). H(k) can, in general have dependence on coordinates 
q . Given the generalized speeds for the joint, the joint map generates the joint linear 
and angular velocity, expressed in the child body frame. The following are used for 
the joints: 



H(k)- 

H{k)- 
H(k)- 

H(k) = 



'A 0 0 0], pin 
0 0 0 A], slider 
~£ 3 Oj], ball 

1 free 



' % 'C*(*)J 

The collection of generalized speeds for all the bodies comprises the 
vector u, the generalized coordinates for the system. As before, access to a function 
that can generate the vector l V k {k) given (q,u) and a specific joint type, is assumed. 
Access to a function that can compute the derivatives q{k) = q{q(k),u(k)) is also 
assumed. This routine generates the time derivative of the generalized position 
coordinates: 

q = W(q)u 

where W(q) is a block diagonal matrix that relates q and u , with each block 
depending upon the joint type: 

q = u for pin j oint, slider j oint 



for ball joint 



where q = \s x e 2 s 4 ] and u = \o\ co 2 co z ] 

and a free joint is a combination of 3 slider joints and one ball joint. Note that there 
are 4 q 's (derivatives of the Euler parameters) associated with 3«'s for ball joints. 

Similarly, l A k (k) , the generalized acceleration of the hinge point of 
body k in its parent, is given by: 







e 2 


_ 1 




2 


A. 





>A k (k)- 



It is these generalized coordinates q, and generalized speeds u, the internal 
coordinates for purposes of this description, of the molecular system which are 
calculated. Rather than working with the typical inertial coordinates (x , y, z) and 
speeds in these inertial coordinate systems, calculations for the subject molecular 
5 system are reduced. 

CALCULATIONS OF THE EQUATIONS OF MOTION 

With the exemplary rigid multibody, torsion angle model described, 
the equations of motion can now be calculated. In accordance with the present 
invention, the motion of the MBS molecular model is determined by the Residual 
1 0 Form. The Residual Form method requires calculations termed the "first" kinematic 
calculations to distinguish them from the "second" kinematic calculations, which are 
further required by the Direct Form (which is included in this description for purposes 
of comparison). 

First Kinematic Calculations for the Molecular Model 

15 In the first kinematic calculations, given the internal coordinates of the 

molecular system, (q,u,u) and the system parameters, the following position, 
velocity and acceleration kinematics are computed for each rigid body k of the 
molecular model. (In passing, it should be noted that when the First Kinematic 
calculations are done for the Residual Form method, the u is passed in as a guess of 

20 the solution which the integration method then refines to the correct solution. In 
contrast, u is set to zero when used for the Direct Form method. This is shown 
clearly in the later descriptions of the two methods.) 
For each body k compute: 

N co\k\ N v Qk {k\V(k), 

V*(*M(Jt) 

25 These computations are done recursively, starting from each base body and 
progressing to the leaves. 

N C k (k) , the direction cosine matrix for body k in ground is defined as: 
*C*(1)='C*(1) 

N C k (k) = N C k {i) i C k {k), k = 2,...n, i = inb(k) 
i C k (k) comes from the joint routine described above. 
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r Q,Qk (k) , the position vector from Q t , the hinge point of the parent of 
body A: to Q k , the hinge point of body k, expressed in the parent frame, is defined as: 

r e - Qk (k) = d(k) + r PkQk (k), k = l,...n 
r PkQk (k) comes from the joint routine. 

r OQk (k) , the position vector from the inertial origin O to Q k , the hinge 
point of body k, expressed in the global frame, is defined 

r OQ k(1) = r Q,Q k(l) 

r OQk (k) = r OQk (i) + N C k (i)r Q& (k), k = 2,...n, i = inb(k) 
l 0 k (k) , the rigid body transformation operator for body k is defined 

V(k) , the spatial velocity for body k at its hinge point, expressed in 
the frame of body k, is defined 

V(k)±(y* {k) } = y k \k)V(i) + k = 2,...n, i = inb(k) 

{ N v & (k)J 

A(k) , the spatial acceleration for body k at its hinge point, expressed 
in the frame of body k, is defined 

where 

^f(M0 + [ ic , w(V(0 X i(0xrgft J 

a?= i C k \k) N a? k (i) 

Of course, the computations can all be computed in a single pass if desired. 

After completing these steps for one incremental time step, the MBS 
can service kinematics requests to compute the (generalized) position, velocity, or 
acceleration information for any point of any body. This is done by computing the 
required information for any point in terms of the hinge quantities for its body, using 
standard rigid body formulas. 
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Residual Computation 

With the first kinematic calculations described above, the residual 
computation for the Residual Form method can be determined. A detailed description 
of the Residual Form and its application to molecular modeling is found in the 

5 previously cited co-pending U.S. Patent Appln. No. , entitled, 

"METHOD FOR RESIDUAL FORM IN MOLECULAR MODELING," filed of even 

date. This computation fills in two partitions of the vector f ^ q given previously. 

\Pu) 

The first partition is called p q , the kinematic residual, and the second partition is 
called p u , the dynamic residual. The kinematic residual is computed from the 
10 difference between a q , which is passed-in from the (implicit) integration submodules 

f!: 66, and the derivative computed by each joint: 

U q-W{q)u = p q 

J* J The dynamics residual is also computed. Starting with a given state of 

W the molecular model, i.e., given (q,u,u) and the system parameters, a program 

p3 15 routine models the 'environment' of the MBS. Such routines are readily available to, 
or can be created by, practitioners in the computer modeling field. The routine takes 
N the values (q, u) determined by and passed in from the integration submodules 66 and 

pj returns (the state-dependent) T(k) = I J > me applied spatial force for a body k at 

its hinge point Q k , and cr(k) , the hinge torque for the body k. The dynamics 
20 residual, p u (k) , associated with generalized speeds u(k) for the body k is then 
computed by the following steps: 

1. Perform the calculations for the molecular model by the Residual Form as 
described above with the passed-in state values (q,u,u) ; 

2. Generate f(k) , the spatial load balance for each body k in the model having n 
25 bodies: 



T(k) = M{k)A{k) H 



N a k (k)( N a k (k)xp(k)) 



-T{k) 



k = \,.. 
3. Compute p u (k) 



-20- 



for A =n to2by-l 

p u (k) = H(k)f(k)-a(k) 

i = inb(k) 

f(i)+=^ k (k)f(k) 
end 

A (l) = tf(l)f(l) 

The Residual Form method evaluates the extent to which the system 
differential equations are satisfied. Zero residual indicates that the applied forces are 
5 in balance with the inertia forces. However, this does not mean the system is in static 
equilibrium, but rather that the applied forces would reproduce the given u when 
applied to the system in the state (q,u) . The residuals can be interpreted as that 
y, additional hinge torque needed to balance the applied and inertia forces. In the 

'i: literature this method is known as either inverse dynamics, or the method of computed 

m 10 torques. It governs the case where the u are all prescribed. At this point all the 
i7i computations required for the Residual Form are complete. The residuals p q and p u 

Z- are used directly by the implicit integrator in the integrator submodule 68. 

Second Kinematics Calculations for the Molecular Model 
t To carry out the Direct Form method, calculations in addition to the 

flj 15 first kinematics calculations are required. These additional calculations are termed 
p the second kinematics calculations. The values P(k),D(k), 'y k (k), 'K k (k) are 

computed as follows: 

1 . Perform the calculations for the Molecular Model by the Residual 
Form as described above, i.e., the first kinematics calculations. 
20 2. P(k) , the articulated body inertia of each body k, is initialized. 

P{k) = M{k), k = l,...,n 
3. The objects below are then generated: 
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for k = n to 2 by -1 

D{k) = H(k)P(k)H\k) 
G = P(k)H\k)D-\k) 
T=£ 6 - GH(k) 

i K k {k) = i 0 k (k)G 
i = inb(k) 

Pd)+=y k (k)P(k)Y\k) 

end 

D(l) = H(l)P(l)H\l) 
The functional dependence of these quantities is only upon the generalized coordinate 
q. Therefore, the first kinematics calculations are programmed in anticipation of 
performing the second kinematics calculations. 

Forward Dynamics Calculations 

Finally, u is calculated by sweeping inboard, then outboard, of the molecule: 
z(k) = 0 6 , k = l,...n 
for& to 2 by -1 

m = p u {k)-H{k)z{k) 

u(k) = D~\k)£{k) 

i = inb(k) 

z(i)+=Y(k)z(k)+ i K k (k)p u (k) 
end 

e(\) = Pu (l)-H(\)z(\) 
u(\) = D-\\)e(\) 
u{l) = v{\) 
S(\)=H\l)v(\) 
for k = 2 to n 
i = inb{k) 

S(k) = y ki, (k)<5(i) + H\kMk) 
u(k) = o(k)- i K k *(k)S(i) 
end 

With the First and Second Kinematics Calculations, and the Forward Dynamics 
Calculations, the Direct Form method is available. 
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Direct Form Method for the Equations of Motions 

The Direct Form method takes the current state (q,u) and computes the 
derivatives (q, u) using the above algorithms, which are then used by the integration 
method to advance time. 
5 Given: (q,u) 

Compute: {q,u) 

1. Compute q using joint specific routine as above 

2. Perform above First Kinematics Calculations with u = 0 

3. Generate residuals p u as above 
10 4. Negate the residuals p u = -p u 

5. Perform Second Kinematics Calculations 

6. Compute u using Forward Dynamics Calculations above 

The Direct Form method produces the hinge accelerations it in 
response to the applied forces acting on the system. Fig. 5 summarizes the 
15 computation steps of the Residual Form method and the Direct Form method. 

JACOBIANS IN IMPLICIT INTEGRATION 

The MD equations which model a molecule (such as a protein), are 
implemented as a multibody system (MBS). These equations represent Newton's 
Laws and are expressed as a set of differential equations y = f{y,i). The differential 

20 equations are implemented using a suite of Order( N ) multibody dynamics methods. 
To advance the equations in time, in accordance with the present invention, an 
implicit method of numerical integration is used, in particular, L-stable implicit 
integration methods, such as implicit Euler, Radau5, and SDIRK3. 

An important ingredient of this integration process is formation of the 

25 Jacobian of the differential equations. This is 

By 

Since the function /is itself computed by an algorithm rather than by an explicit 
formula, the Jacobian computation represents a substantial amount of work. In the 
simplest approach, the Jacobian can be formed numerically by differencing the 
30 derivative routine. This is a delicate operation because the quality of the Jacobian is a 
tradeoff between round-off and truncation errors. Typically half the working 
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precision in the result is retained by choosing a good perturbation size in the 
difference scheme. In practice, though, this is difficult to do. 

However, the structure of the governing equations may be exploited to 
improve the Jacobian computation. The exemplary multibody dynamics methods 
5 illustrate this. The algorithms involved compute exact derivatives, even though 

numerical methods are used to execute the formulas. The derivatives obtained are in 
error by amounts that depend upon round off and the conditioning of the multibody 
system under consideration. But no approximations are involved at the equation 
level. 

10 In general, G, the iteration matrix used in the Newton loop of the 

implicit integrator has the form G = E-aJ , where E is the identity matrix and a is 
be some scalar function of the timestep. See the previously referenced U.S. Patent 

p Appln. No. , entitled "METHOD FOR LARGE TIMESTEPS IN 

O MOLECULAR MODELING," filed of even date, for a description of implicit 

7. J 15 integrators. Changes in step size require refactoring G , but not reforming J . 
W Reforming / is needed only when the Jacobian is needed at a new state. G is used 

p in a linear subproblem within a Newton loop. The following is solved: 

U GAy = -r(y i n ), 

where r(y) is the residual function for that particular implicit integration method, 
p 20 As shown later, J has a special structure, which is inherited by G . 

This means that equation solving with G can be done at reduced cost, if the structure 
is exploited. 

The quality of the Jacobian affects the ability to solve the nonlinear 
equations resulting from discretization in the integrator. Failure to solve the Newton 
25 loop may require retraction of a trail step and reduction of the integration time step. 
The timestep should be controlled by accuracy, rather than failures in the Newton 
loop. 

Computation of Analytic Jacobian 

The Jacobian J is a matrix which represents a linearization of the 
30 equations of motion. Normally, the governing equations for a dynamical system are 
linearized around an equilibrium state, or perhaps a state of steady motion. In this 
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J{q,u) = 



case, the equations are linearized around an arbitrary state so all possible contributing 
terms should be developed. It is customary to describe Jin terms of its partitions: 

'dq d£ 
dq du J J gg J qi 
du du J m 
K dq du ; 

Structure of J qq AND J qu 

The q equation is q - Wu , where the matrix Whas a block-diagonal 
structure. Each block depends upon the joint type. Pins and slider joints give rise to 1 
by 1 identity blocks. A ball joint generates a 4 by 3 block that expresses the Euler 
parameter derivatives in terms of Euler parameters and body angular velocity measure 
numbers (generalized speeds). 

From the q equation above, these two partitions of the Jacobian matrix 

are formed: 

These equations are to be interpreted for symbolic purposes only. In 
practice, there is no need to generate the matrix W explicitly. The non-zero block 
diagonal elements are filled in as discussed in the previous section on the kinematic 
residual. 

Structure of J m AND J uu 

The u derivatives are more complicated. Since u = -M~ 1 p u , 
d(-M~ x p u \ d(-M-'p\ 

— '- and — ^ L must be computed. (Note that p u is the residual of the 

dq du 

dynamics routine developed earlier). Here a key result from the field of Numerical 
Analysis is used to avoid the derivative of the matrix inverse. 

Suppose A(x)y(x) = b(x) , then we can write y{x) = A~ l {x)b{x) . If 

y(x 0 ) = z is known, and the value of ^ X ^ \ x=x0 must be obtained, we have 



dy\ 
cbcL 



= A- 1 (xO) — (b(x)-A(x)z)\ 

) dx L 
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where z = A~ x b is held fixed when forming the right-hand side of the equation. 
Applied to the described multibody routines, 

§1 = -±(M->p u (q,u,Q)) = -M~ l ■Z-p m iq,u,z) 

dx ox v ' ox 

where z = M~ l p u (q,u,G) . This result avoids the computing of the derivative of M' x , 
which is a hypermatrix. The matrix inverse is "pulled-out". In the above equation, 
"x" is either the generalized coordinate q or the generalized speed w, depending on the 
Jacobian partition that is to be computed. 

In summary, to compute the blocks J uq and J uu , three steps are 

followed: 

1 . Given (q, u) compute z = -M' l p u (q, u, 0) using the Direct Form method. 
This simply says to compute u from the current state and save it as the 
variable z. 

2. Compute the analytic Jacobian of the dynamics residual routine. In this step, 

the matrix —p u (q,u,z) is to be formed. This quantity clearly depends upon 

dx 

the vector z computed in Step 1. Note that the numerical value of the residual 
p u in this step is zero for each element since we are computing the Jacobian 
around a consistent solution of the motion equations. The partitions J uq and 
J uu of the Residual Jacobian are obtained by substituting "q" and "u" for the 
"x" above. 

3 . Back-solve the result of Step 2 with the mass-matrix to obtain the desired 
block. The back-solve operation is accomplished in the Direct Method routine 
by processing a residual vector into a u vector. The Second Kinematics Step 
only needs to be performed once, since the back-solves are done at the 
nominal value of the state. In fact, the Second Kinematics routine must have 
been called in Step 1 while computing z, so the variables should still be 
cached. 

In words, the Jacobian of our derivative routine can be formed by back-solving the 
Jacobian of our residual routine. The residual Jacobian is much simpler to compute 
than the derivative Jacobian. Steps 2 and 3 above are derived in the following 
sections. 
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The Residual Jacobian 

The computation of the Residual Jacobian is closely related to the 
Residual Form method for dynamics, which is summarized here: 

1 . Perform an outboard pass that computes the kinematic data that depends upon 
5 position and velocity. 

2. Make a call to the force routine which generates the atomic forces and 
consolidates them into spatial loads acting on the bodies. 

3. Perform another kinematic pass that computes acceleration level quantities 
(using passed-in u ), and combines inertia forces with the spatial loads from 

10 Step 2. 

4. Perform an inward pass that generates the residual at each joint. This pass 
recursively computes the resultant of the (spatial) inertial and applied forces 
for the 'nest' of bodies outboard of the joint in question. The residual is the 
projection of the resultant on the joint's degrees of freedom, given by the joint 

15 map data. 

At a high level the residual computation can be considered to depend 
upon two kinds of forces: 'motion forces' and external forces. The motion forces are 
computed directly by the multibody system. The external forces are available to the 
multibody system from a force modeling routine that computes the various 

20 interatomic forces such as electrostatics and solvents. A similar procedure is followed 
when computing the Jacobian. The multibody system builds the Jacobian of the 
motion forces, and combine it with the Jacobian of the external forces. 

Partitioning the Force Field into Effects 

There are many forces that may be acting on the molecule, and these 

25 forces may be computed in various intrinsic coordinate frames that are most 

convenient for that particular force "effect". For example, electrostatic terms may be 
computed using multipole methods and spherical coordinates, covalent terms may be 
computed in terms of torsion and bond angles, and solvent forces may be computed in 
global Cartesian coordinates. During the computation of the Residuals, these forces 

30 are transformed from their intrinsic coordinate frame to the MBS coordinates. 

The same exchange occurs to compute Jacobians. The native 
Jacobians in their intrinsic coordinates are be brought into the MBS coordinates. This 
requires the use of the chain-rule to transform between intrinsic and the MBS 
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generalized coordinates. It is important that each effect co-computes its function 
value and Jacobian, because many of the same terms are needed for each 
computation. Each effect is transformed into a set of spatial loads T effect {k) , where k 
is the index of a generic body in the system. The totality of these effects is given the 
5 symbol T(k) . 

Effect Jacobians Brought into the Residual Jacobian 

At a high level, the residual routine was previously implemented from 

the equation 

p u =H0(MA-T) 

1 0 The implementation uses Order( N ) methods which are immediately obvious from the 
equation above. In this equation T is a vector of spatial loads acting on the pivots of 
f a the multibody system, where each element is a spatial load (a 6-vector composed of 

p. one force and one torque). It actually represents all effects other than inertia loads or 

U-* pure hinge loads. The term in parentheses represents the load balance for each body, 

jjj 1 5 The first term is the inertia force, the next is the spatial load. M(k)A(k) is the spatial 

inertia force for a typical body. This is built from the body mass properties and the 
f spatial acceleration of the body pivot. The spatial acceleration is computed before the 

U residual routine is executed by the Forward dynamics routine. The operator H <f> is 

i implemented in a routine that performs an Order( N ) inward pass. 

H 20 Even without knowing anything about the details of the computation 

implied by the equation above, the contribution of — , the spatial load Jacobian (the 

effect Jacobian) to , the residual Jacobian, can be immediately inferred: 

dx 

SA.-H^+... 

dx dx 

The ... refers to terms not involving the effect Jacobian. Again, q or u for "x", is 
25 substituted, depending upon which partition of the Jacobian is being computed. 

The role of an effect Jacobian in the residual Jacobian is the same as 
the role of the effect in the multibody equations. This means that T contributes to the 

residual p in the same way that a column of — contributes to . Both are 

ox ox 

processed by the same operator H0 . This is a crucial point because it means that no 
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new method is needed for this part of the Jacobian computation. (A different method 

i • dT. 
is need to obtain — — ). 

dx 

Thus, given the effect Jacobian, its contribution is assembled into the 
residual Jacobian by operating on it with the original residual routine, treating it as a 
5 multi-column set of spatial load vectors. This is a direct consequence of the linearity 
of the equations. 

The columns of the residual Jacobian play the same role in the 
derivative Jacobian routine as the residual vectors play in the Forward Dynamics 
routine. The dynamics routine performs a back-solve on a data vector it receives, and 
1 0 doesn't need to know what the data is, just what operation to perform on it. This 
applies to all the routines. 

h* • Adding the . . . terms, the chain rule is used to show the whole 

O 

p equation: 

ti 



d Pu _ ( d(MA) dT ^ | g 8W , d(Hy) 
dx I dx dx J dx dx 



1 5 At this level, there are four contributions to the Jacobian: the inertia forces, the spatial 
forces, and contributions due to changes in the operators 0 and H. The quantities z 
and y refer to (MA - T) , and <f>z , respectively, which are held fixed while evaluating 
the last two terms. The numerical values of these terms are already available from the 
residual computation. Another observation about the above equation is that the 

20 operator ^ depends only upon q, but not u. Thus, this term remains constant while 

computing the partition J uu . Similarly, the spatial load can be split into its constituent 
effects, some of which do not depend upon u. In general, this means that knowledge 
of the multibody equations can be exploited to optimize the computations. 

Up to now, what to do with — once the term has been computed has 

dx 

25 been described, but a description of how to form the term has not been made. These 
details are in the following sections. 

Computing the Effect Jacobians and Combining with the Residual Jacobian 

So far a high-level description of the Jacobian computation has been 
given. It can be seen that the computation has a very algorithmic flavor to it. There 
30 are very distinct phases to the task, just as there were also for the Forward Dynamics 
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routine described previously. There, computation of atomic forces is clearly the 
bottleneck step. Yet the overhead in the multibody equations for dealing with forces 
is fairly small. In that case, a call was made to the force routine, and what occurs 
inside the routine was ignored. When it comes to the Jacobian, this aspect is less true. 

A call is still made to obtain the Effect Jacobian, but there is a lot of 
processing needed before the Effect Jacobian can be assembled into the Residual 
Jacobian. The details of dealing with force fields to produce Jacobians are covered in 
the next sections and an example of incorporating electrostatics is developed. All 
other loads follow a similar development. 

Electrostatics as an Example Effect 

The basic premise of electrostatics is that the force between two 
charged particles is 

p - qiq j - 

This is a classical inverse square law. F tj , the force acting on particle i due to particle 
j, depends upon the charge of the particles, attractive for oppositely charged particles, 
repulsive for like charges. The symbol r & is a unit vector directed from particle i to 
particle j; r y is the distance between the particles; a: is a unit-dependent constant 
related to the strength of the forces. Of course, the force acting on particle j is equal 
and opposite to that acting on the particle i. Hence, given a collection of particles 
interacting through Coulomb's Law given above, the net force acting on each particle 
is computed by summing the pair-wise forces. For this example, the forces are 
computed in global Cartesian coordinates. 

With the atomic forces in hand, the multibody forces can be generated. 
The system of forces acting on the particles of each body is replaced by a spatial load 
acting at the pivot of each body. The atomic forces are first expressed in a body-fixed 
basis, and then shifted to the pivot using the station coordinates of the particular atom 
to which the force is bound. 

Fy is a vector. The derivative of F & with respect to dr y , small changes 
in the particle's relative positions is D y , a tensor. It is such that dF i} = Dy»dr & . For 
Coulomb force the tensor is 
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Some observations: 

1 . The Force Jacobian is a matrix of size n atoms x n atgms . Each element is a 3 by 3 
tensor. The (ij) block gives the derivative offeree on atom i with respect to 
small changes in the position of atom j. In general, every force model is 
required to support an intrinsic Jacobian method for analytical processing. 

2. Storage requirements for the Force Jacobian quickly become impractical. This 
leads to the notion of interface "contraction" where the Jacobian of all the 
forces acting pair-wise between the atoms are reduced or "contracted" to 
acting pair-wise between the bodies. 

3. Because the Coulomb force is a pair- wise interaction, each force contributes to 
two blocks in the overall Jacobian. Thus, each force is processed at constant 
cost, and the overall Jacobian is computed at a cost proportional to the number 
of atoms squared, i.e., Order( TV 2 ). This is the same as the computational cost 
of the force itself! This is a rather good result for computing analytic 
Jacobians. A numerical Jacobian requires a fresh force computation each time 
an element of the state is perturbed. This leads to cubic growth, i.e., 

Order( iV 3 ), in the cost of the numerical Jacobian. Hence, the analytic 
Jacobian is much cheaper to compute as well as more accurate than a 
numerical Jacobian. 

4. The computation of the Jacobian is conveniently done in the same routine that 
computes the force (co-computation). However, it typically needs to be done 
far less often than the force computation. Therefore, a flag can be used to 
trigger the Jacobian calculation only when needed. 

Coupling to the Displacement Gradient 



process the Jacobian further. This is due to the fact that the multibody system is 
formulated using relative coordinates. The chain rule is applied to each atomic force 
F(k,i) and is called "coupling to the displacement gradient". This denotes the global 
Cartesian force on atom i, which resides on body k. 



Having obtained the (intrinsic) Force Jacobian, it is necessary to 



8F(k,i) 
dqj 




dF(k,i) dr(p,s) 
8r(p,s) d qj 



p=\ szp 
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where r(p, s) is the position of atom "s" on body "p". 

The first term in the sum selects an element of the Force Jacobian 

which was just computed. The quantity dr( f ,s ^ is an element of the displacement 

oq j 

gradient. A typical term gives the change in an atom's position due to a small change 
in a generalized coordinate. Note that this term is strictly a kinematical quantity 
having nothing to do with the force computation. Thus, the Force Jacobian can be 
computed once and then continually reprocessed by the chain rule for each coordinate 

drj . 

in the multibody system. This step represents a matrix vector multiply, since — is a 

column vector with n atom entries (each a 3 vector), and the Jacobian is a square 
matrix n atoms x n atoms , where each element is a 3 by 3 tensor. 

It is possible to improve this computation, since many of the entries in 
the displacement Jacobian are known to be zero. This is due to the fact that 
incrementing a particular hinge does not displace every atom in the system, but only 
those outboard of the displaced hinge, and not disjoint from the hinge in question. 
For instance, rotating the base body induces a change in all atoms of the system. But 
perturbing a torsion angle on any terminal body induces a change only to those atoms 
resident in the terminal body. Therefore, roughly speaking, about half the work may 
be saved by optimizing the computation. This reduction comes from a strictly 
knowledge-based approach. 
Interface Contraction 

In the process of forming T(k), the spatial load on body k, the load comes from the 
atomic forces acting on atoms that belong to body k. Each force is transformed from 
global to local coordinates, and then shifted to the body pivot. A concise statement of 
this procedure is: 



T(k) = ^(k,i)T(k,i) 



iek 



The operator <fi(k,i) is 
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where p{k,i) are the fixed station coordinates of atom i on body k. Note that the new 
quantity T(k,i) appears. This is just the atomic force turned into a spatial load at the 
atomic position of atom i: 



The first element of the atomic spatial load is zero because there is no torque exerted 
by the force field on individual atoms. 

Now, T(k) relates atomic forces to body spatial loads. So, the 
derivative of this equation relates differential atomic forces to differential spatial 
loads: 

dqj tT dqj d qj 

= T t (k) + T 2 (k) 

The second term, T 2 (k) , in this equation is discussed at the end of this section, as it 
involves the spatial loads, but not the load derivatives. This means the term can be 
treated generically, without worrying how the spatial loads were computed. 
Substituting the definition of T(k), into T x (k) : 

dT(k) dr(p,s) 
~Uhdr{p,s) dqj 

where now the symbol dT<<k) = T>(£,0 dT ^ l \ is an element of the row-reduced 
dr(p,s) tk dr(p,s) 

Force Jacobian. This Jacobian relates differential (body) spatial loads directly to 

differential atomic displacements. Again, r(p,s) refers to the global position of the 

atom s on the body p. The term is formed by summing each atomic Force 

dr(p,s) 

Jacobian element into the destination element in the reduced Jacobian, weighted by 
the atoms' 0{k,i) matrix. Each element of the row-reduced Jacobian is a 6 by 3 
matrix. Hence the rows of the Force Jacobian have been contracted. The contraction 
is evident in the notation: the numerator has only a body index, while the denominator 
has both a body and an atom index. Depending on the number of atoms per body, the 
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row reduction can provide a savings in both storage and execution time when 
differential spatial loads must be formed. 

Note that the row-reduction procedure only needs to be done once 
before computing the Residual Jacobian. The overhead of performing the reduction is 
more than offset by the reduced cost of the smaller matrix vector products which must 
be formed. 

8T(k) 



Note that in forming 



dr{p,s) 



there is no need to save the elements of the atomic 



Force Jacobian. That is, each element d3"(ft»0 only nee( j s to ^ e ava ilable while its 
dr{p,s) 

contribution to ^1^1 i s being computed. So, more than one element of the big 
dr(p,s) 

Force Jacobian is not required at a time. 

The global coordinates of a typical atom, r(p,s) , are computed in 
terms of r(p) , the global coordinates of body p's pivot, and p(p,s) , the station 
coordinates of the atom: 

r(p,s) = r(p) + N C k (p)p(p,s) 
By differentiating we find, (with some results from the kinematics of finite rotations): 

^^ = ^^ C \p)p(p,s))A(p) 
d qj dqj 

Augmenting this equation with the additional equation A(p,s) = A(p) and defining w, 
the spatial derivative 



20 the result is: 





' A~ 




w = 


dr 
Pi. 






' A(p,s) ' 






w(p,s) = 


dr(p,s) 


= <? 




dr{p) 




L 




d qj 



The vector 1 can be interpreted as generating a rate of change of orientation for each 
body. It is a field quantity, in the sense that it can potentially vary at each point in 
25 space. For rigid bodies undergoing pure rotations (without deformation), it is constant 
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The second term, T 2 (k) , is given by: 



and involves the original spatial loads T(k) and derivative of the operator 0(k,i) : 

VN C k (k) 0, j 
O3 N C k (k) 



5 Thus 



\E 3 p(k,i)l 



fe* |_ Ms 



w rfC*(*) p\k,i) N dC k (k) 



N dC k (k) 



T(k,i) 



N dC k (k) = N dC k (i) 'C*(Ar) + *C*(z) 'dC*(A) 
10 is computed recursively from the base body outward and 

>dC k (k) = ?f^-dq(k) k = \,...n 
dq{k) 

where dq(k) is defined in the next section, and - - — , the partial derivative of the 
interbody direction cosine matrix is a function of the joint type connecting the bodies: 



pin : 



rc k (k) 



= - Js 3 sin(q k ) + A cos(q k ) + XX sin(^ ) 



slider : — = 0, 



ballandfree: 8 ' C ^ k K 2 



d l C k (k) 



0 


s 2 






-2s x 


-s 4 


<?3 




-2s, 



8^C k ik}_ 



d l C k (k) _ 
ds 4 



--2s 2 


3 


s 4 


s x 


0 


s 3 


-s 4 




-2s 2 


--2s, 


-s 4 


£ x 


s A 


-2s 3 




s \ 


s 2 


0 


' 0 




s 2 ' 


s 3 


0 




-s 2 


s \ 


0 _ 
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In the next Section the Force Jacobians are combined with the Inertia Force Jacobians to 
finally form the Jacobian of the Residual Routine. 

The Residual Jacobian 

dT(k) 

The previous Section describes the formation of the Force Jacobian , 

dw(p) 

5 which must be coupled to the spatial displacement gradient, in order to form the derivative of 
the spatial forces. This section describes the formation of the spatial displacement gradient 
and the formation of the Jacobian of the residual routine. 

The recursive algorithms for computing the entire Jacobian are described. The 
Jacobian algorithm is not actually set up to compute the Jacobian. As is typical of automatic 
10 differentiation routines, it computes the matrix vector product J uq dq + J uu du for arbitrary 
y, passed-in values of the vectors dq and du. In practice, to compute the Jacobian, the "Jacobian 
5 Routine" is effectively called repeatedly with a series of Boolean vectors (a vector with one 
Ift entry set to 1 and all other entries set to zero.) Each call generates the corresponding column 
f% of the Jacobian. Note that some of the steps have already been or are being computed for the 
4-1 5 Residual Form method or the Direct Form method (the Forward Dynamics Calculations), but 
l" are reproduced here for clarity. 

fT 1. Given (q,u) compute z = -M~ ] p u (q,u,0) using the Direct Form method. Also set 

O 

r '~ u = z and recompute A{k) , then p u , which recomputes T(k) . 

~rl 2. Perform contraction to compute the fully row- and column-reduced Force Jacobian, 

™ d^(^) as d escr ibed in section, "Interface Contraction": 



dw(p) 



dT ^dT 
dw dr 



Steps 3 through 10 below are used to fill the columns of J uq : 

3. Compute the analytic Jacobian partitions of the q terms: 

dq 

using joint routines similar to those needed for the First Kinematics Calculations. 

4. Compute q derivatives of position quantities and terms for spatial gradient: 
Previously described methods were used to fill in certain joint-specific fields. These 
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quantities consisted of i C k {k) , the interbody direction cosine matrices, r Q,Qk (k) , the 
spanning vector for each body, and H(k) , the joint map for each body's inboard joint. 
To refer to the derivative of each of these quantities, a prefix 'd' is added to the 
symbol name to make this reference generically. Thus, l dC k {k) means the derivative 
of the direction cosine matrix l C k {k) . 

Each interbody direction cosine matrix (and all the joint-specific) quantities depend 
only on the generalized coordinates of an individual joint. Thus, l dC k (k) is nonzero 
only when the derivative is taken with respect to any of the coordinates for body k. 
To properly 'seed' the recursions being studying, a vector dq is passed in to the 
routine. For Jacobian computation we set one entry is set to 1, and all the other 
entries to 0. Then, the needed preliminary quantities are generated by a typical loop: 

i dC k (k) = ^^-dg(k) k = l,...n 
The partial derivatives of the direction cosine matrices are generated analytically and 
displayed in the section, "Interface Contraction" above. These partial derivatives do 
not depend upon the particular column of the Jacobian that is being computed. 
Setting a particular entry of dq to 1, and all the rest to 0 has the effect of annihilating 
the correct subset of the seed quantities. 
8r &Qk (k) 



d( lk 



-, the partial derivative of the interbody spanning vector is given by 



dr aa (fr) 
dr e ' & (k) 
dr aQ *(k) 

dr e& (k) 



= A{k), slider 
= 0 3x4 , ball 
= 0 a 



pin 








0 0 


0 


1 0 


0" 


0 0 


0 


0 1 


0 


0 0 


0 


0 0 


1 



free 



A(k) here refers to the body's sliding axis which connects it to its parent body. 



dHjk) 
8<2k 



the partial derivative of the joint map is 
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With the above definitions of the partial derivatives, the recursions are seeded with 
the following loop: 

for k = 1 to nbod 
i dC^k) = °^-dq(k) 

dr^(k)= 8rQ ' Qk(k) dm 

end 

After execution of these loops, all bodies have l dC k (k) , dr aQk (k) , and dH(k) , their 
interbody derivative quantities available. 

One new quantity needed in the spatial displacement gradient computation is 
also computed. This is A(k) from the section on Interface Contraction, the rotation 
axis that generates the rate of change of orientation for each body outboard of a 
perturbed joint. Here, this variable is given the symbol d6(k~) , the differential 
rotation axis for each body is 

for k = 1 to nbod 

d0(k) = A{k)dq{k), pin 

d0(k) = %, slider 

d0{k) = not needed for ball, free 

end 

Since arbitrary perturbations to a set of Euler parameters do not produce a pure 
rotation, column contraction cannot be used when computing the corresponding 
column of the Jacobian. The row -reduced Force Jacobian can still (and must) be 
used.. 

After seeding the recursions, N dC k (k) , dr 0Qk (k) , { d0 k (k) , dA(k) are computed: 
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N dC k {\)= i dC k {\) 
dr° & (l) = dr e ' B *(\) 
l d<f> k {\)=% 
dJi(X) = d0{\) 
for k = 2 to nbod 

N dc k (k) = N dc k (iyc k (k) + N c k (iydc k (k) 

dr 0Qk (k) = dr 0& (0 + N dC k (i)r Q& (k) + N C k (i)dr e& (k) 



dZ{k) = 'C'WdZii) + dd{k) 

end 

where 

a = l dC k {k), b = dr e& (kyC k (k) + r && l dC k {k\ 
c = %,d = i dC k {k~) 

Compute q derivatives of velocity: 

A loop that computes the rate of change of joint velocity due to change in joint angle 
starts the process: 

for k = 1 to nbod 

l dv k (k) = dH\k)u(k) 

end 

This quantity is the rate of change of joint velocity due to change in joint angle. 
Obviously, it is nonzero only for joints whose map contains coordinate dependence. 
For free joints, the generalized speeds produce relative linear velocity that depends 
upon the joint orientation. 

After computing ! dv k (k) , dV(k) , the derivative of the spatial velocity of each 
body, is computed. This is done by the following loop: 

dV(V)= i dv k {\) 
for k = 2 to nbod 

dV(k) = W/V(0+ i 0 k *dV(i)+ i dv k (k) 
end 

Couple Force Jacobian to spatial displacement gradient to compute 7] (k) 
for k = 1 to nbod 

"^ dT(k)dw(p) 
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Compute second term of the Force Jacobian T 2 (k) and append to T t (k) : 
fork= 1 to nbod 

~ N dC k (k) p(k,i) N dC h (kj\ 
N dC k (k) 



dT{k) = r,(*) + £ 



T(k,i) 



end 



Compute q derivatives of acceleration-related terms: 

Again the process starts with a loop that computes l da k {k) = dH*{k)u(k) : 

for k = 1 to nbod 

i da k (k) = dH\k)u(k) 

end 

This is the change in joint acceleration due to a change in coordinate. Then, dA(k) , 
the derivative of the spatial acceleration of each body is computed. 
dA(l) = 'da k (l) 

_'v Q *(k)\ l'dv 3 *(k) 



where -> means the 6 vector is decomposed into two 3 vectors, 
for k = 2 to nbod 
dtl = l dC k * (k) ( N co k (i) x ( N a> k (i) x r && (k))) + 
\ N d<v k (i)x( N co k (i)xr Q - & (k)j) + 
l C k * (*) ( N w k (/) x( N da> k 0) x fOA (Jt) )) 4 
( *V (/) x ( N co k (0 x dr Q& (*))) 

da = '<//*(*).4(0 + ! /*(jt)dL4(0 + ^J 
c0j± i C k \k) N a k (i) 

dtl = dcOjx l CO k (k) + Q}jx l dcO k (k) 

dt3 = 2<to y x 'v& (A:) + 2^ x Wv a (£) 

r^2i 



[^3 



1 • * 
\+'da*(k) 

J 
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The symbols introduced here with = are meant to be temporary variables not needed 
after computation of dA(k) . 

After computing the spatial acceleration derivatives, the computation of 
df{k) , the spatial inertia force derivatives, is performed: 

for k = 1 to nbod 

dvl = N dco k (k) x Jg t (k) • N G) k + N (O k (k) x Jg t (k) • N dco k 

dv2 = N dco k (k) x ( N co k (k) x p(k)) + N a> k (k) x ( N dco k (k) x p(k)) 

df(k) = M(k)dA(k) + j - dT(k) 

end 



9. Compute d p(k) , the joint residual derivative for body k: 

for k = nbod to 1 

dp(k) = dH(k)f(k) + H(k)df(k) - da{k) 
i = inb(k) 
if i >0 

df(i) = df(i) + wwnk) + '/(k)df(k) 

end 
end 

After executing this routine, the values stored in dp(k) are the new column of the 

Residual Jacobian — p u (q,u,z) . 
dq 

10. Back-solve the — result of previous step with the mass-matrix to obtain the desired 

dq 

dq' 




The back-solve operation is accomplished in the Direct Form method routine by 
processing a residual vector into a u vector. The Second Kinematics Calculations 
only needs to be performed once for the whole Jacobian, since the back-solves are 
done at the nominal value of the state. In fact, the Second Kinematics routine must 
have been called in Step 1 while computing z, so the variables should still be cached. 
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Steps 1 1 through 13 below are used to fill the columns of J uu : 

11. Compute u derivatives of velocity: 

This routine takes a passed-in vector du and computes 'dv k (k) = H*(k)du(k) . Then, 
dV(k) , the derivative of the spatial velocity of each body, is computed: 

dVQ) = 'dv k (l) 

for & = 2 to nbod 

dV(k) = ^ k \k)dV(i) + 'dv k (k) 

end 

12. Compute the velocity-induced derivative df(k) . As presented here, the routine is 
specialized for the case of no velocity dependent external loads. The surviving terms 
are those due to changes in inertia forces alone. Even if there were changes in 
external loads, it would only be required to include the contribution of dT(Jc) as 
before. 



»-[t] 



dA(l) 

for k = 2 to nbod 

[(W(/)x(V(/)x^(i))) 

(W(0x(W(0xr oa (/t)))^ 



dt\= l C k \k) 



da± i 0 k \k)dA(i) + ^ i 

coj^'C^ikYo/ii) 
dcoj=' l C k \kyd<Q k {i) 
dt2 = da>jx i co k (k) + a>jx l d co k (k) 
dt3 = Ida j x ' v a (k) + 2a?jx l dv Qk (k) 

end 

After computing the spatial acceleration derivatives, df(k) , the spatial inertia force 
derivatives, is computed: 
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for k = 1 to nbod 

dvl = N dco k (k)^l Qk (k) ■ N o) k + N a> k (k) x ^ (k) • N da? 
dv2 ± N dco k (k) x( N a k (k)x p(k)) + N co k {k)x( N dm k (k) x p(*)) 

df(k) = M(k)dA(k) + 
end 

13. Compute d p(k) , the joint residual derivative for body k: 

for k = nbod to 1 
dp(k) = H(k)df(k) - dcr(k) 
i = inb(k) 
if i >0 

df(i) = df(i)+ i /(k)df(k) 

end 
end 

After executing this routine the values stored in dp{k) are the new column of the 

Residual Jacobian — p u (q,u, z) . 
du 

14. Back-solve the — result of previous step with the mass-matrix to obtain the desired 

du 

8u_ 
du " 

du du 

The back-solve operation is accomplished in the Direct Method routine by processing 
a residual vector into a u vector. The Second Kinematics Calculations only need to 
be performed once, since the back-solves are done at the nominal value of the state. 
In fact, the Second Kinematics routine must have been called in Step 1 while 
computing z, so the variables should still be cached. 

The above steps complete the computation of the analytic Jacobian as long as 
the forces only have dependence on q . This accommodates the classical situation where all 
atomic forces are derivable from a potential function. To accommodate velocity-dependent 
forces, such as simple viscous damping, some of the above steps need to be modified as 
follows: 
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In Step 2 above, we also need to compute the contracted velocity Jacobian 

d-^) wn i c h is block diagonal, must also be computed. 
dr(k) 

In Step 6 above, the computation of T x (k) must be augmented with the 
contracted velocity Jacobian: 

for k = 1 to nbod 

f dT(k) dwjp) | y 8T(k) dr(k,i) 
]{ } Udw(p) dqj £rar(M) dqj 
end 

where 

= N dc k (*) [ N v & (*) + W (*) x 0] + 

S ,V C* (it) [ A V/v a (it) + W (A) x 0] 

Ifi 

yj A step is then added after Step 11, which is called Step 11a. This new step 

-110 computes dT(k) : 

f* t?dr(k,i) duj 

rj where 

S f^M= jv c*(it)rw t (it)+ w ^(ik)x/<it,oi 

5m, L J 

While executing Step 12 above, the last loop for df(k) is modified by 
1 5 subtracting the velocity-dependent force derivative dT(k) : 
for & = 1 to nbod 

dvl = N daf (k) x ^ (k) ■ N co k + N eo k (k) x ^ (k) ■ N daf 
dv2= N d& k (k)x( N a k (k)xp(k)) + N a k (k)x("da k (k)xp(k)) 

dT(k) = M(k)dA(k) + 11- dT{k) 
end 

The rest of the Steps remain the same. 

Fig. 6 summarizes the operational steps of the Analytic Jacobian method, 
which has been described in detail above. 
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Fig. 7 shows a plot of the accuracy of the numerical Jacobian versus the 
accuracy of the analytic Jacobian for an exemplary MD system. In the best case in which the 
perturbation was perfectly selected, the digits of accuracy for the generalized coordinates (q) 
and generalized speeds (u) from the numerical Jacobian, illustrated by line 152, were still 
5 only half that of the analytic Jacobian, illustrated by line 1 50. 

ADDITIONAL EMBODIMENTS 

The present invention has many embodiments besides the examples described 

above. The list below has other embodiments and applications: 

• Order of Forces included in Jacobian 
10 Any order of the forces to be included in the Jacobian, include, but not limited 

to Order( N ), Order( N 2 ), Order( N 3 ), and Order( TV 4 ). An example of an Order( N ) force 
L ; field would be an electrostatic force field using fast multi-pole expansion methods (see, for 
P example, Greengaard, The Rapid Evaluation of Potential Fields in Particle Systems, 
If! Massachusetts Institute of Technology Dissertation, 1 988) rather than direct method which is 
jjjl5 Qrder(JV 2 ). 

J;; • Analytic Jacobian for Direct Form 

a When the governing equations are in Direct Form, the so-called "forward 

l\ dynamics" form of the equations is obtained. In this form, the equations process a state 
O vector and applied efforts and generate the acceleration at each of the joints modeled in the 
p20 system. 

* u = M~ l (/) 

The Jacobian then represents the partial derivatives of the accelerations with respect to 
elements of the state vector. The preferred embodiment shows several algorithmic methods 
for computation of these partial derivatives. The methods are exact and do not utilize 
25 numerical approximations to form derivatives. 

The Direct Form produces the u partitions of the Jacobian: 
S(M-'(f)) 
dq 

and 

a(M-(/)) 

du 

30 by using an algorithmic counterpart to the function which computes the u function. 
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